function [SE_alpha, SE_variety] = std_error(alpha_hat, variety_hat, ...
    moment_mean_all, var_cov, W, N_star,...
    K, D, N,...
    pr_fc, x_fc, lapse, beta_c, crra, rho, sigma_c1, sigma_c2, eulerc,...
    N_y, pr_y, y_medi, T1, T2, B1_c, B2_c, naic_unique,...
    no_insurer_fe, ind_select, delta, use_log_n, myopic,...
    resource_y1, resource_y2)
%% preliminaries
% number of moments used
N_moments = length(moment_mean_all);

% number of parameters estimated internally
num_est = length(ind_select);

% order of x = [alpha; lapse_u; variety_u; xi];
x = [alpha_hat; 0; variety_hat; zeros(N.i,1)];

% store sample moments evaluated at the estimates
moment_mean_bench = moment_mean_all;

%% calculate the derivative of moments
moment_deriv = zeros(N_moments, num_est); 

scale=1e-7;
for k=1:num_est
    
    if k==1 % alpha
        deriv = abs(alpha_hat)*scale;
        x_deriv = [alpha_hat+deriv; 0; variety_hat; zeros(N.i,1)];
    elseif k==2 % variety_u
        deriv = abs(variety_hat)*scale;
        x_deriv = [alpha_hat; 0; variety_hat+deriv; zeros(N.i,1)];
    end
    
    [~, ~, ~, ~, ~, ~, ~, ~,...
        ~,~,...
        moment_mean_all_delta, ~,...
        ~, ~] ...
        = demand_objective_inner(x_deriv, K, D, N,...
        pr_fc, x_fc, lapse, beta_c, crra, rho, sigma_c1, sigma_c2, eulerc,...
        N_y, pr_y, y_medi, T1, T2, B1_c, B2_c, naic_unique,...
        no_insurer_fe, ind_select, delta, use_log_n, myopic,...
        resource_y1, resource_y2);

    moment_deriv(:,k) = (moment_mean_all_delta-moment_mean_bench)/deriv;
end

%% calculate the asymptotic variance
temp1 = inv(moment_deriv'*W*moment_deriv);
V = temp1*(moment_deriv'*W*var_cov*W'*moment_deriv)*temp1;

SE_alpha   = sqrt((1/N_star)* V(1,1)); 
SE_variety = [];

% display
disp(' ')
disp('****************')
disp('DEMAND ESTIMATION RESULTS')
disp('Parameter   Estimate    Std     Estimate/Std')
fprintf('%-12s %10.6f %10.6f %10.2f \n', 'alpha', alpha_hat, SE_alpha, alpha_hat/SE_alpha);


